A minimal integer automaton behind crystal plasticity 
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Power law fluctuations and scale free spatial patterns are known to characterize steady state plastic flow in 
crystalline materials. In this Letter we study the emergence of correlations in a simple Frenkel-Kontorova (FK) 
type model of 2D plasticity which is largely free of arbitrariness, amenable to analytical study and is capable 
of generating critical exponents matching experiments. Our main observation concerns the possibility to reduce 
continuum plasticity to an integer valued automaton revealing inherent discreteness of the plastic flow. 



At the macroscale one usually assumes that crystalline ma- 
terials flow plastically when averaged stresses exceed yield 
thresholds. At the microscale plasticity evolves through a se- 
quence of slow-fast events involving collective pinning and 
depinning of dislocational structures. Classical engineering 
theory has been very successful in reproducing the most im- 
portant plasticity phenomenology such as yield, hardening 
and shakedown [IJ, however, a fully quantitative link between 
the phenomenological theory and the microscopic picture of 
plasticity remains elusive. The main reason is that the phe- 
nomenological approach implies spatial and temporal averag- 
ing in the system with poorly understood long range correla- 
tions. 

The presence of such correlations have been confirmed 
by numerous experiments revealing intermittent character of 
plastic activity with power law statistics of avalanches and self 
similar structure of dislocation cell structures |2|. The emer- 
gence of power laws suggests that in plasticity the relation 
between the microscopic and the macroscopic models is more 
akin to turbulence than to elasticity |[3]|4]|. Similar critical fea- 
tures of stationary nonequilibrium states have been observed 
in a variety other driven systems with threshold nonlinearity 
and rate independent dissipation including, for instance, tec- 
tonic faults, magnets, and superconductors |5 |. The origin of 
scale free attractors in such systems is a subject of active re- 
search, in particular, the problem of classifying the universal- 
ity classes remains largely open |6|. In this situation finding 
the minimal representation of each class which is amenable to 
rigorous analysis is of significant general interest. 

The experimental evidence of plastic criticality has been 
corroborated by several numerical models |7|. The two main 
microscopic approaches are discrete dislocation dynamics, ac- 
counting for dislocation interactions on different slip planes 
ISHTOl, and a pinning-depinning model dealing with plastic- 
ity on a single slip plane UJj l12J. Different meso-scopic 
continuum models implying partial averaging have also been 
shown to generate power law statistics of avalanches with re- 
alistic exponents |4, 13 1. Since scale free dislocation activity 
is expected to be independent of either microscopic or macro- 
scopic details, one can try to maximally simplify the underly- 
ing physics while still capturing the observed exponents and 
even characteristic shape functions |14|. Presently the only 
analytically tractable models of plasticity are the mean field 
theory 1 15 , 16 1, the FRG models of elastic depinning 1 17 1 and 
the Abelian automata of sand pile type 1 1 8 1 . In this group only 
the automata models have a potential of capturing the whole 



complexity of the dislocation patterning and the goal of this 
Letter is to propose a formal reduction of a realistic plasticity 
model to a spin model with discrete time evolution. Instead 
of straightforward time discretization of continuum dynamics 
1 19] we search for inherent temporal discreteness hidden be- 
hind the conventional gradient flow dynamics 1201 . 

While the ID automata, describing successfully plastic hys- 
teresis and rate independent dissipation fall short of capturing 
plastic criticality |20|, the 2D automaton based models may 
be already adequate at least for FCC and hexagonal crystals; 
it is also noteworthy that intermittency has been mostly ob- 
served under single slip conditions 1 15 , 21 1. In such cases one 
can get a realistic model by assuming that plasticity proceeds 
through the motion of a set of parallel edge dislocations. We 
further neglect vectorial nature of the problem and reduce the 
crystal to an array of coupled FK chains |22|. In contrast to 
more conventional DDD modeling | 8 1 where nucleation and 
propagation rules are not associated with the same thresholds, 
the FK type models describes adequately both the multiplica- 
tion of the defects and their kinetics including the finite size 
of the Peierls stress 1231 . 

To describe the model we first recall that classical contin- 
uum dislocation mechanics deals with the energy ^{u) = 
/ 4>(Vu) where u{x) is the displacement field and the func- 
tion 4> is quadratic. The displacement field is allowed to have 
finite discontinuities [u] whose evolution is governed by phe- 
nomenological kinetic relations 1241 . The atomic structure 
of dislocations can be addressed by introducing an internal 
length scale a (Burgers parameter) and replacing continuum 
energy by the discrete one which can be schematically rep- 
resented as ^{u) = a ^([ixj/a) where now u is 3. lattice 
field, [u] is a discrete increment and the function (j) is peri- 
odic with infinite symmetry group |25|. Assuming that the 
order parameter is scalar, one can minimize out remaining 
linear strain variables: the simplest example of the resulting 
dressed description is the one dimensional FK model |26|. 
Our minimal 2D setting can be viewed as an array of cou- 
pled FK chains with the energy ESlEl ^{u) = Y^ij <t>{0, 0, 
where 0{i^ k) = u{i + 1, /c) — u{i^ k) is an axial strain and 
^(i, k) = u{i^ k -\- 1) — u{i, k) is a shear strain; the displace- 
ment field u{i^k) is defined on slNxN square lattice (a = 1). 
The potential (j) is assumed to be quadratic in and periodic 
in ^ (see Fig. [T]); to avoid synchronization we also added 
quenched disorder ^(i9,0 = g{0 + fi^f - ^i^ - ^20. 
where /ii,2(^, j), are independent gaussian random variables. 
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Figure 1 : Periodic dependence of energy on shear strain. Each min- 
imum corresponds to a shp between adjacent atomic planes by one 
lattice spacing. Dislocations correspond to the boundaries of slipped 
regions. Red: continuous potential. Blue: Piece-wise quadratic po- 
tential with the same linear elastic moduli. One elastic domain is 
shadowed. 
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Figure 3: Probability distribution of dissipated energy in a steady 
state; an insert shows the structure of fluctuations during a typical 
cycle. 



Observe that the variables and ^ are not independent and the 
long-range interactions can be revealed through minimizing 
out the non order parameter variable ll36ll . 

The dynamic equations are taken in the form of the over- 
damped gradient flow uii = —d^{u)/du^ where u is the ratio 
of the internal time scale and the time scale of the driving 
11201 . The driving in shear is performed through displacement 
controlled boundary condition Xl^o^ (.{h^) = t, where t is 
the slow time playing the role of loading parameter; in lon- 
gitudinal direction we assume periodic boundary condition 
SiLo^ 6>(i, k) = 0. In our numerical experiments we took 
i^^ = 2, TV = 512 and g{^) = (27r)-2(l - cos(27rO). The 
initial state was dislocation free and the dispersion of disorder 
varied in the range 0.01 — 0.2. For computations we used an 
implicit-explicit FFT method 1271 . 
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Figure 2: (a) Macroscopic strain-stress curve showing plastic hys- 
teresis and shakedown. An insert illustrates a fragment of the de- 
formed lattice with two dislocation dipoles nucleated around an im- 
perfection: red and blue colors correspond to dislocations of differ- 
ent sign, (b) Evolution of dislocation density with cycling. An insert 
shows small scale fluctuations. 

The results of direct numerical simulations (automaton re- 



duction is discussed later) are presented in Fig. |2(a)| where we 
show the macroscopic strain- stress curves covering first two 
cycles of loading/unloading in the hard device. Notice that 
hysteresis loops converge indicating that the system exhibits 
plastic shakedown. Reaching steady state is marked by the 
stabilization of dislocation density, which also shows a char- 
acteristic nucleation related overshoot (see Fig. 2(b)| ). Steady 
state yielding is characterized by the formation of stable dis- 



location structures (cells) with plastic activity limited to in- 
termittent dislocational exchanges between the clusters; the 
latter remain mostly stable from one cycle to another but have 
a finite lifetime as in observations |28 1. To separate individual 
avalanches we introduce an irrelevant threshold and define the 
avalanche energy by integrating viscous dissipation over its 
duration: E = N~'^ ^ • j J v?dt. We observe that the prob- 
ability distribution P{E) stabilizes after several cycles (see 
Fig. |3]) exhibiting a robust power law behavior with expo- 
nent e ^ 1.6 ± 0.05 obtained by maximum likelihood method 
1291 . This value is in perfect agreement with experiments in 
ice crystals and fits the generally accepted range 1.4-1.6 |2l|9l; 
most remarkably it is also consistent with the value obtained 
for 2D colloidal crystals 1301 . The approximate proportion- 
ality between the plastic slip size and the dissipated energy 
ensures that acoustic emission measurements would exhibit 
the same exponent e. 

The spatial counterpart of the observed time correlations 
is the fractal structure of dislocational patterns. The disloca- 
tion rich regions (clusters) can be identified by the localized 



peaks of the energy density landscape (see insert in Fig 4(a) ) 
and the corresponding probability density shows a power law 
structure with exponent 1.45 ±0.1 (Fig |4(a)| ). Another way 
to quantify the fractal clustering is to compute the correlation 
function of the dislocation distribution C{r) ~ ISTl . We 
observe that during the first loading cycle D ~ 2.0, which 
is expected given the random nature of the quenched disor- 
der. With cycling the long range correlations develop (see Fig. 
4(b)|and in the shakedown regime we record D ~ 1.74 inde- 



pendently of initial disorder. Note that dislocation patterns 
with I) ^ 1.64 — 1.79 have been observed experimentally in 
crystals with multiple slip systems; in simulations with a sin- 
gle slip system fractal patterning has been previously linked 
to the possibility of dislocation multiplication [|32i which is 
operative in our model. 

Despite the conceptual transparency of the above model, 
the mechanism of reaching the critical regime remains ob- 
scure. The model, however, can be simplified further if we 
notice that in quasi-static limit ^ the relaxation is instan- 
taneous and the system remains almost always in equilibrium 
d^/du = 0. In order to solve equilibrium equations analyti- 
cally we can replace the smooth periodic potential by a piece- 
wise quadratic potential (see Fig. [T]) defined in each period 
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Figure 4: (a) Probability distribution of dislocation rich regions in 
the shakedown state; an insert shows spatial distribution of the en- 
ergy density (j){0^^)\ (b) Correlation function C(r) after the first 
cycle (squares) and after the fifth cycle (circles); an insert shows a 
characteristic stress field during steady yielding. 

{{d - 1)^^ {d + 1)^0) as g{£,) = ^(^ - where d is the 
new integer- valued spin variable describing a quantized slip; 
since at given lattice field j) the equilibrium equations are 
linear the strain field can be found analytically (cf. 1 35 1). The 
Fourier image ^(q) of the shear strain reads 

|(q) = (4(q)s,-(q)d(q) + JT(q))/A(q), 

where q = (qx^Qy) = {2iTk/ N ^2itI / N) is the wave 
number. Here we defined i^(q) = 5~(q)5+(q)/ii(q) + 
5^(q)4(q)^2(q) and A(q) = 2i^(cos(g,) - 1) + 
5-(q)s+(q) where 5j(q) = ±(1 - cos((7a) ± ^ sin((7a)) with 
a = x^y. Notice also that we control the average shear strain 
^^o(q) = ^^(q) and that the quenched disorder becomes the 
source of the residual strain ^/^(q) = i^(q)/A(q). 

It is now straightforward to reformulate the model as an 
integer valued automaton. Observe that the variable = 
^ ~ (^0 + representing shear strain fluctuations must be 
confined between the thresholds —(^ — ^(/i, t) < A^(i, j) < 

- ^{h, t), where ^(/i, t) = [^/,]-^ + t and [']-^ denotes the 
inverse Fourier transform. When reaches the threshold, 
the integer parameter d is updated d ^ d -\- M(A^), where 

r +1, ifA^>e-ah,t), 
M{Ao = l -1, if < -e - ah,t) 

[ otherwise. 

After each increment of loading t we recheck the stability un- 
til all the units are stabilized; the dissipated energy during an 
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Figure 5: Scaling collapse of the dissipated energy in automaton 
model for different system sizes N. An insert displays the kernel 
L{x,y) and illustrates the update process in the automaton. 

avalanche is the difference of the total energies for two subse- 
quent stable states. 

The use of cellular automaton representation greatly re- 
duces the complexity of numerical computations while the 
behavior of the system remains the same including the shape 
of the stress strain hysteresis, the evolution of the dislocation 
density and the structure of spatial and temporal correlations. 
To illustrate the statistics we show in Fig. [5] the finite-size 
scaling collapse of energy dissipation at the critical state; here 
we assumed that P{E) = E~^(p{E/Ec) with universal cut 
off function (p and the cut off energy which diverges in the 
thermodynamic limit as Ec ^ . Our computations show 
that again e ^ 1.6 ± 0.05 and that 6 ^ 1.2 ± 0.1 which 
is close to the value ^ ^ 1 obtained for plastic strain incre- 
ments in (gl [33| . These exponents are insensitive to the de- 
gree of disorder in the studied range; for larger disorder we 
observed a cut-off which is no longer size dependent. Based 
on our computations we conclude that both models, the one 
with continuous dynamics and smooth potential, and the one 
with discrete dynamics and piece wise quadratic potential be- 
long to the same universality class. This behavior is markedly 
different from the prediction of the mean field theory where 
smooth and cusped potentials lead to different universality 
classes |34|. 

An important question is whether the toppling rules in our 
automaton are Abelian meaning that the outcome of the in- 
stability in multiple sites does not depend on the toppling or- 
der. The update of the "slope" field A^ can be represented in 
Fourier space as A^ A^ — I/(q)M(A^), where 



sm{qy/2Y^Ksm{q,/2Y 

is the analog of the toppling matrix in the sand-pile mod- 
els. The corresponding kernel L{x^y) in the real space is 
highly anisotropic, long-range and conservative (see insert in 
Fig. [5]). We compared numerically all conventional updating 
strategies and found that the microscopic configuration shows 
some small dependence on the choice of the strategy while the 
macroscopic observables, including the shakedown hysteresis 
loop and the statistics of avalanches (critical exponents) re- 
main unaffected. One can conclude that our automaton has a 
weak (statistical) form of Abelian symmetry which may still 
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be helpful for the mathematical analysis fTS]. Another im- 
portant property of our automaton model is that it necessar- 
ily lowers the energy during each avalanche. The dissipative 
structure is obvious in the continuum model and is inherited 
by the automaton model. 

In conclusion, to elucidate the origin of self organized crit- 
icality in plasticity we reduced a realistic continuous dynam- 
ics to an integer automaton by replacing the fast dissipative 
stages with jump discontinuities controlled by random thresh- 



olds. The fact that despite the long range character of elas- 
tic interaction the computed exponents are different from the 
predictions of the mean field theory may mean that at least for 
some crystal classes plasticity is effectively a 2D phenomenon 
laying below the upper critical dimension. 
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